Error Diffusion Using Linear Pixel Shuffling

Linear pixel shuffling error diffusion is a digital halftoning algorithm that combines the linear pixel shuffling (LPS) order of visiting pixels in an image with diffusion of quantization errors in all directions. LPS uses a simple linear rule to produce a pixel ordering giving a smooth, uniform probing of the image. This paper elucidates that algorithm. Like the Floyd-Steinberg algorithm, LPS error diffusion enhances edges and retains high-frequency image information. LPS error diffusion avoids some of the artifacts (“worms,” “tears,” and “checkerboarding”) often associated with the Floyd-Steinberg algorithm. LPS error diffusion requires the entire image be available in memory; the Floyd-Steinberg algorithm requires storage proportional only to a single scan line. 1. Error Diffusion Halftone Production The error diffusion algorithm transforms a gray scale image, I, with pixel values in the interval [0:0; 1:0℄, to a black-andwhite image, B, with values in f0; 1g. The following pseudocode describes error diffusion: for every pixel position i,j in I if I[i][j] < 0.5 then B[i][j] = 0 else B[i][j] = 1 error = I[i][j] B[i][j] distribute the error among unprocessed neighbors of i,j The order of pixel visitation generally takes the form of raster processing: for ( i = 0; i < i_max; i++ ) { for ( j = 0; j < j_max; j++ ) { process pixel i,j } } The input and output images dimensions are i max by j max. Floyd and Steinberg’s error diffusion algorithm [4] follows this pixel ordering and distributes the error to four unprocessed neighbors of I[i][j] according to the kernel DFS : DFS = 1 16 P 7 3 5 1 (1) Figure 1: Floyd-Steinberg error diffusion. Notice the “worm” artifacts at the ends, and “tearing” in the middle. (All the ramp images are 64 256 pixels.) where P denotes the pixel currently being processed. I[i][j+1] += error * 7/16 I[i+1][j-1] += error * 3/16 I[i+1][j] += error * 5/16 I[i+1][j+1] += error * 1/16 This process produces bilevel images with visual appearance capturing the full range and detail of the original image. This is particularly effective in case the original image has a lot of detail. The resulting images do, however, often contain “worm” artifacts in very dark and very light regions, and a “tearing” or “checkerboarding” where the image’s original gray value was slowly varying around 0:25, 0:5, or 0:75; see Fig. 1 (our ramp figures are presented here using enlarged pixels–approximately 80dpi—in order to illustrate the results of the algorithms and avoid transformations caused by unknown printing processes; the proper viewing distance is 4–5 feet). Floyd-Steinberg error diffusion was also used to render the “Woman with French Horn” in Fig. 2 (approx. 130dpi). 2. Linear Pixel Shuffling Linear pixel shuffling (LPS) is a method for visiting image pixels distributed evenly all over an image. LPS uses a simple linear rule, given by matrix multiplication (modulo a parameter depending on the image size) i j = M xy (2) The original image processing pseudocode uses Eq. (2), becoming: Figure 2: Floyd-Steinberg error diffusion of “Woman with French Horn” (image courtesy of Heidelberger Druckmaschinen A. G.). for ( x = 0; x < N; x++ ) { for ( y = 0; y < N; y++ ) { multiply M * (x,y)’ to obtain (i,j)’; process pixel i,j } } The matrix M and the parameter N are described below. Define the Fibonacci-like sequence G by the recurrence: G0 = 0 G1 = 1 G2 = 1 Gn+1 = Gn +Gn 2 for n 2 (3) Terms G0 through G14 of this sequence are 0; 1; 1; 1; 2; 3; 4; 6; 9; 13; 19; 28; 41; 60; 88 We also need this sequence with negative subscripts (the definition allows us to work backwards in the subscripts easily). Terms G 1 through G 14 are 0; 1; 0; 1; 1; 1; 2; 0; 3; 2; 3; 5; 1; 8 For our algorithm, we round the image size up to a square of side Gn (this is the N in the above pseudocode) and ignore pixels that fall outside the actual image. Define a Gn Gn table T with entries defined as Tpq = (pGn 2 + qGn 1) (mod Gn) (4) Because g d(Gn 2; Gn 1; Gn) = 1 (5) for all n, every number in 0; 1; 2; ; Gn 1 occurs in T exactly Gn times. Here is a portion of the 88 88 table for 0 p; q < 13. 0 60 32 4 64 36 8 68 40 12 72 44 16 41 13 73 45 17 77 49 21 81 53 25 85 57 82 54 26 86 58 30 2 62 34 6 66 38 10 35 7 67 39 11 71 43 15 75 47 19 79 51 76 48 20 80 52 24 84 56 28 0 60 32 4 29 1 61 33 5 65 37 9 69 41 13 73 45 70 42 14 74 46 18 78 50 22 82 54 26 86 23 83 55 27 87 59 31 3 63 35 7 67 39 64 36 8 68 40 12 72 44 16 76 48 20 80 17 77 49 21 81 53 25 85 57 29 1 61 33 58 30 2 62 34 6 66 38 10 70 42 14 74 11 71 43 15 75 47 19 79 51 23 83 55 27 52 24 84 56 28 0 60 32 4 64 36 8 68 A particularly nice feature of the table T is that values which are numerically close are physically distant. Examine a 7 7 block centered at one of the zeros in T : 49 21 81 53 25 85 57 2 62 34 6 66 38 10 43 15 75 47 19 79 51 84 56 28 0 60 32 4 37 9 69 41 13 73 45 78 50 22 82 54 26 86 31 3 63 35 7 67 39 Notice that the smaller numbers are not close to 0. Because T was constructed using a simple linear rule, this phenomenon holds for any value, not just zero; it is even more visible for larger values of Gn. The LPS pixel selection algorithm works by first visiting the elements in the Gn Gn square that correspond to the positions (i; j) for which Tij = 0, followed by positions (i; j) for which Tij = 1, and so on. We may express this in the processing pseudocode: for ( x = 0; x < N; x++ ) { process all N pixels i,j such that T[i][j] = x } The matrix M that achieves this is M = G n+1 Gn 3 G n Gn 2 (6) For example, using N = G14 = 88 M = G 13 G11 G 14 G12 = 1 28 8 41 (7) 2.1. Mathematical Justification of the LPS Algorithm The above algorithm processes all Gn Gn image pixels for the following reason. Let i j = G n+1 Gn 3 G n Gn 2 xy (8) We will show that Tij = x (9) Multiplying out the above matrix formula, we have Tij = iGn 2 + jGn 1 (10) = (G n+1x+Gn 3y)Gn 2 +(G nx+Gn 2y)Gn 1 (11) = (G n+1Gn 2 +G nGn 1)x +(Gn 3Gn 2 +Gn 2Gn 1)y (12) = (G n+1Gn 2 +G nGn 1)x (13) +(Gn 3 +Gn 1)Gn 2y (14) The values in Eq. (8–14) are all modulo Gn. One of the factors in (14) satisfies Gn 3 +Gn 1 = Gn 0 (mod Gn) (15) by (3); so Tij is independent of y. The coefficient of x in (13) satisfies Gn 2G n+1 +Gn 1G n 1 (mod Gn) (16) for all n, because Gn 2G n+1 +Gn 1G n +GnG n+2 = 1 (17) for all n; this is easily established by mathematical induction. Finally, because three consecutive terms of the G sequence are relatively prime (5), M 0 y takes Gn distinct values for 0 y < Gn. Linearity allows us to conclude that i j ranges over all Gn Gn pixel positions. 2.2. Other Applications of LPS There are a wide variety of applications for LPS; most of them are explained in detail in [2]. Here are several that pertain directly to image processing. LPS was originally developed to aid development of computer graphics image synthesis, especially fractal creation (see [1]). By rendering pixels in the LPS order, a low resolution indication of the eventual picture is almost immediately visible. Slow computers or ambitious programmers otherwise produce images that can take hours to see—starting with the extremely boring topmost rows of pixels. Figure 3: An 88 88 image with the first 968 (12.5%) of the pixels marked black. The pixel ordering for image rendering suggests an image file format. The initial portion of such a file gives a low resolution image; subsequent portions fill in the details. Because LPS is linear, the pixels that are rendered earliest can be magnified (“fat pixels”) to cover the entire image area very early. Subsequent pixels can be magnified by smaller factors, eventually not magnified at all. If there are constant patches in the image, the later pixels often make no difference to the image appearance—they do not need to be rendered or stored at all. This observation leads to a lossless image compression technique for text and line-drawing images that is competitive with run-length encoding but has the useful property that a prefix of the compressed file is a lossy compression of the entire image. The results of image morphology operations (dilation, erosion, opening, closing) can be rapidly approximated by performing the basic operations in LPS order. For some applications, an approximation that takes, say, 15% of the total time may be sufficient. The table T can be used as a halftone mask. Fig. 3 shows an 88 88 image with the first 968 of the pixels marked. This is the dot pattern associated with an LPS mask at gray level 0.125 (where 1.0 indicates black). Images can be searched to locate objects (“Where’s Waldo?”) or gather image statistics, such as histogram estimation. An infinite version of LPS is appropriate for Monte Carlo integration. Figure 4: The ramp rendered using LPS error diffusion. 3. LPS Error Diffusion By visiting the pixels uniformly all over the image, LPS allows us to diffuse the errors in all directions. For example, we can use the kernel which John Szybist investigated for his Computer Science Master’s Project [3]: DLPS = 1 32 266664 1 1 1 1 2 3 2 1 1 3 P 3 1 1 2 3 2 1 1 1 1 377775 (18) and distribute the error of the processed pixel P to 20 of its neighbors. In order not to loose any of the dispersed error, we need to keep track of any neighborhood of a zero value in T (a 5 5 neighborhood in this case) and note the positions in the kernel corresponding to already quantized pixels. We thus update the diffusion kernel after each completed inner loop (the yloop in our pseudocode), potentially replacing some non-zero coefficient with a zero and updating the divisor (initially 32 here) correspondingly. In addition to the Szybist filter (18) used for Figs. 4 and 5 we can use a wide variety of kernels; several are shown in Fig. 6. Ramps rendered using these kernels are shown in Fig. 7. 4. Comparisons and Conclusions LPS error diffusion does eliminate worms, tearing, and checkerboarding artifacts. Unfortunately, LPS halftoned images appear to be mottled in comparison to Floyd-Steinberg halftoned images. The LPS halftoning algorithm requires that the entire image be present in memory, in order to diffuse error in all directions. Floyd-Steinberg only requires one or two scan lines (depending on coding cleverness, to diffuse errors only to the right and down). The number of operations in both cases is linear in kernel size and number of image pixels. The texture variations of images rendered by various LPS halftoning kernels may present opportunities for achieving good dot patterns of ink on paper that eludes most current digital halftone systems. Figure 5: Linear pixel shuffling error diffusion.